#ifndef __TDP2_RANDOM_HPP__
#define __TDP2_RANDOM_HPP__

#include "tdp2.hpp"

namespace TDP
{
	namespace Math
	{
		template<class _Tfloat,LargeInt _Dim> class Vec;
	};

	class Randomizer
	{
	private:
		UInt32 _state[624];
		UInt32 _index;
	public:
		inline void Seed(UInt32 seed)
		{
			_state[0] = seed;
			_index = 624;
			for(unsigned int i=1;i<624;i++)
				_state[i] = UInt32(0x6c078965 * (_state[i-1] ^ (_state[i-1] >> 30)) + i);
		}

		inline Randomizer(UInt32 seed)
		{
			Seed(seed);
		}

		inline Randomizer()
		{
			Seed(0x5088ff07);	// pi * 10^10 with first nibble dropped
		}

		inline void RebuildState()
		{
			for(int i=0;i<624;i++)
			{
				UInt32 y = (_state[i] & 0x80000000) | (_state[(i+1) % 624] & 0x7fffffff);
				_state[i] = _state[(i+297) % 624] ^ (y >> 1) ^ ((y&1)*0x9908b0df);
			}
			_index = 0;
		}

		inline UInt32 RandomUInt32()
		{
			if(_index == 624)
				RebuildState();

			UInt32 y = _state[_index++];

			// Temper
			y ^= y >> 11;
			y ^= ((y << 7) & 0x9d2c5680);
			y ^= ((y << 15) & 0xefc60000);
			y ^= y >> 18;

			return y;
		}

		inline Float32 RandomFloatNot1()
		{
			return static_cast<Float32>(RandomUInt32() & 0x7fffff) / (8388608.0f);
		}

		inline Float32 RandomFloat()
		{
			return static_cast<Float32>(RandomUInt32() % 0x800001) / (8388608.0f);
		}

		inline Float64 RandomDouble()
		{
			Float64 baseline = (Float64)(RandomUInt32() & 0x3ffffff) * 67108864.0;
			baseline += (Float64)(RandomUInt32() & 0x3ffffff);
			return baseline / 4503599627370496.0;
		}

		inline bool RandomBool()
		{
			return (RandomUInt32() & 1) != 0;
		}

		class Math::Vec<Float32,3> RandomDirection();
		class Math::Vec<Float32,3> RandomDirectionLambertian(const class Math::Vec<Float32,3> &normal);
		class Math::Vec<Float32,3> RandomPointOnTriangle(const class Math::Vec<Float32,3> &p1, const class Math::Vec<Float32,3> &p2, const class Math::Vec<Float32,3> &p3);
		class Math::Vec<Float32,2> RandomPointInCircle(const class Math::Vec<Float32,2> &center, Float32 radius);

		class Math::Vec<Float32,3> RandomDirectionHigh(Float32 distance = 1.0f);
	};
};
 

#include "tdp2_math.hpp"
#include "tdp2_math_angle.hpp"
#include <math.h>

namespace TDP
{
	inline Math::Vec<Float32,3> Randomizer::RandomDirection()
	{
		Float32 u = RandomFloat() * 2.0f - 1.0f;
		Float32 n = sqrtf(1.0f - u*u);

		Math::Angle<Float32> theta = Math::Angle<Float32>(2.0f * Float32(Math::Pi_k) * RandomFloatNot1());

		return Math::FVec3(n * theta.Cos(), n * theta.Sin(), u);
	}

	inline Math::FVec3 Randomizer::RandomDirectionHigh(Float32 distance)
	{
		Float64 d = RandomDouble() * 2.0 * Math::Pi_k;
		Float64 s = sin(d) * Float64(distance);
		Float64 c = cos(d) * Float64(distance);

		Float64 u = RandomDouble() * 2.0 - 1.0;
		Float64 n = sqrt(1.0 - u*u);

		return Math::FVec3(Float32(n * c), Float32(n * s), Float32(Float64(distance) * u));
	}

	inline Math::FVec3 Randomizer::RandomDirectionLambertian(const Math::Vec<Float32,3> &normal)
	{
		Math::FVec3 side(normal[2], -normal[0], normal[1]);
		side = (side - normal *side.DotProduct(normal)).Normalize2();
		Math::FVec3 side2 = side.Cross(normal);

		Float32 rf = RandomFloat();
		Float32 u = sqrtf(rf);		// Square root = rate of photons at N = proportional to N, Lambertian
		Float32 n = sqrtf(1.0f - rf);

		Math::Angle<Float32> theta(2.0f * Float32(Math::Pi_k) * RandomFloatNot1());

		return side * (n * theta.Cos()) +
			side2 * (n * theta.Sin()) +
			normal * u;
	}

	inline Math::FVec3 Randomizer::RandomPointOnTriangle(const Math::FVec3 &p1, const Math::FVec3 &p2, const Math::FVec3 &p3)
	{
		Float32 v1 = sqrtf(RandomFloat());
		Float32 v2 = RandomFloat() * v1;

		return p1 + (p2 - p1)*v1 + (p3 - p2)*v2;
	}

	inline Math::FVec2 Randomizer::RandomPointInCircle(const Math::FVec2 &center, float radius)
	{
		Float32 dist = sqrtf(RandomFloat()) * radius;
		Math::Angle<Float32> angle(2.0f * Float32(Math::Pi_k) * RandomFloatNot1());

		return Math::FVec2(center[0] + angle.Cos()*dist, center[1] + angle.Sin()*dist);
	}
}

#endif
